function K = stiff_bb_2(dofs,V,T,d,Int_UU_cell,Int_UV_cell,Int_VV_cell)
%function K = stiff_bb_2(dofs,V,T,d,Int_UU_cell,Int_UV_cell,Int_VV_cell)
% This is hp version, where d is a vector
%
%comput the total striff matrix
if nargin == 4    
    d_max = max(d);
    Int_UU_cell = cell(d_max,1);
    Int_UV_cell = cell(d_max,1);
    Int_VV_cell = cell(d_max,1);
    for degree = 2:d_max
        Mat = build(degree-1,degree-1);
        Id = diag(ones((degree+1)*(degree+2)/2,1));
        Du = dirder(Id,1,0,-1); % the direction derivate on v2-v1
        Dv = dirder(Id,0,1,-1);  % the direction derivate on v3-v1
        Int_UU_cell{degree} = Du'*Mat*Du;
        Int_UV_cell{degree} = Du'*Mat*Dv;
        Int_VV_cell{degree} = Dv'*Mat*Dv;
    end
end

%% calculate areas for each triangle
x21 = V(T(:,2),1)-V(T(:,1),1);
y21 = V(T(:,2),2)-V(T(:,1),2);
x31 = V(T(:,3),1)-V(T(:,1),1);
y31 = V(T(:,3),2)-V(T(:,1),2);
areas = 0.5*(x21.*y31 - x31.*y21);

% finish parameter pareparing
m = 0; nt = size(T,1);
for k=1:nt
    degree = d(k);
    md = (degree+1)*(degree+2)/2;
    m = m + md^2;
end
Indx1 = zeros(m,1);
Indx2 = zeros(m,1);
S = zeros(m,1);
pos_start = 1;
for k = 1:nt
    degree = d(k);
    Int_UU = Int_UU_cell{degree};
    Int_UV = Int_UV_cell{degree};
    Int_VV = Int_VV_cell{degree};
    V1=V(T(k,1),:);V2=V(T(k,2),:);V3=V(T(k,3),:);
    %area = abs(triarea(V1,V2,V3));
    x13 = V1(1)-V3(1);
    y13 = V1(2)-V3(2);
    x23 = V2(1)-V3(1);
    y23 = V2(2)-V3(2);
    LocK = ((x23^2 + y23^2)*Int_UU - ...
            (x13*x23 + y13*y23)*(Int_UV+Int_UV') +...
            (x13^2 + y13^2) * Int_VV )/abs(4*areas(k));   % J = 2*tri_area(K)
    [i,j,s] = find(LocK);
    loc_dof = dofs(k,1):dofs(k,2);
    L = length(i);
    Indx1(pos_start:(pos_start + L-1)) = loc_dof(i);
    Indx2(pos_start:(pos_start + L-1)) = loc_dof(j);
    S(pos_start:(pos_start + L-1)) = s;
    pos_start = pos_start + L;   
end
dim = max(max(dofs));
K = sparse(Indx1(1:(pos_start-1)),Indx2(1:(pos_start-1)),S(1:(pos_start-1)),dim,dim);
